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We extend to nonequilibrium processes our recent theory for the long time dynamics of flex- 
. ible chain molecules. While the previous theory describes the equilibrium motions for any bond 

■ or interatomic separation in (bio)polymers by time correlation functions, the present extension of 

pL^ ' the theory enables the prediction of the nonequilibrium relaxation that occurs in processes, such as 

, T-jump experiments, where there are sudden transitions between, for example, different equilibrium 

1 states. As a test of the theory, we consider the "unfolding" of pentadecane when it is transported 

from a constrained all-trans conformation to a random-coil state at thermal equilibrium. The time 
evolution of the mean-square end-to-end distance (-RendM)noneq after release of the constraint is 
computed both from the theory and from Brownian dynamics (BD) simulations. The lack of time 
translational symmetry for nonequilibrium processes requires that the BD simulations of the relax- 
^ ■ ation of (-R^nd W)noneq be computed from an average over a huge number of independent trajectories, 

rather than over successive configurations from a single trajectory, which may be used to generate 
equilibrium time correlation functions. Adequate convergence ensues for the non-equilibrium sim- 
ulations only after averaging 9000 trajectories, each of 0.8 ns duration. In contrast, the theory 
' requires only equilibrium averages for the initial and final states, which may be readily obtained 

from a few Brownian dynamics trajectories. Therefore, the new method produces enormous savings 
in computer time. Moreover, since both theory and simulations use identical potentials and solvent 
O | models the theory contains no adjustable parameters. The predictions of the theory for the relax- 

ation of (i?cnd W)noneq agree very well with the BD simulations. This work is a starting point for 
I ' the application of the new method to nonequilibrium processes with biological importance such as 

, the helix-coil transition and protein folding. 

t> ■ 
r- : 

CN . I. INTRODUCTION 

<N ; 

qq ■ The nonequilibrium dynamics of polymers is important for understanding a variety of phenomena. For instance, 
the dynamics of protein folding/unfolding is a fundamental biological process. Our understanding of this process , 
however, is severely limited by a shortage of adequate theoretical methods for studying the protein dynamics with 
realistic molecular models. Experimentally, it is difficult to achieve sufficiently short time resolution to follow the 
full dynamics of protein folding/unfolding. Conversely, the main theoretical technique, molecular dynamics, cannot 
reach time scales long enough to probe many of the interesting events due to the tremendous amount of computation 
required. 

We have been developing ^ theory for the long time dynamics of flexible polymers and polypeptides_Ln solution. 
Our mode coupling theoryErB is based on a significant extension of the generalized Rouse-Zimm theoryH~£3 to include 
contributions to the dynamics from the memory function ( "internal friction" ) terms that are customarily ignored. Prior 
papers demonstrate that the mode coupling theory produces an excellent representation of the long time equilibrium 
time correlation functions for both alkanes and short polypeptides. The present paper describes the extension of 
the mode coupling theory to treat the nonequilibrium dynamics occurring when a system in equilibrium (or in a 
constrained equilibrium) is suddenly transported to a different equilibrium (perhaps constrained) state, such as in 
T-jump experiments. 

We apply the theory to describe the relaxation of pentadecane from nonequilibrium initial conditions to equilibrium 
at a temperature of 300 K. The pentadecane chain begins in its maximally stretched conformation, in which all the 
dihedral angles are trans. The molecule subsequently relaxes to equilibrium, where there is a distribution of dihedral 
angles due to thermal fluctuations. The direct simulation of this non-equilibrium process requires orders of magnitude 
more computer time than is necessary for simulating equilibrium time correlation functions. The complexity of 
the non-equilibrium simulations is the reason pentadecane is taken to begin in the all-trans configuration, rather 
than using an initial thermal distribution at a temperature lower than 300 K, where a huge number of additional 
non-equilibrium trajectories would be necessary to average over the thermal distribution of initial conditions. In 
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contrast, the treatment of such a T-jump process by the theory merely doubles the more limited computational effort 
necessary to describe the dynamics beginning from the all-trans pentadecane. Studying the "unfolding" dynamics 
of this relatively simple system provides insight into how to apply the mode coupling theory to more complicated 
biological systems, such as the unwinding of an a-helix or the unfolding of a protein. 

A crucial step in the theory is the selection of a suitable basis set. Due to symmetry, the basis functions for describing 
nonequilibrium dynamics must differ from those previously used for calculating equilibrium dipole time correlation 
functions. However, the same principles apply independent of the dynamical observables of interest. For equilibrium 
dipole time correlation functions, a simple first-order basis set produces reasonable estimatesEI □ of alkane dynamics, 
andpajiiore sophisticated sorting method yields excellent approximations while maintaining a manageable basis set 
sizeJj'tj The present paper extends our mode coupling procedure to determine the appropriate basis sets for describing 
the long time portions of the nonequilibrium dynamics. As in previous studies, we test the effectiveness of various 
basis sets by comparing the predictions of the theory with "exact" results from Brownian dynamics simulations. 
Because all the equilibrium inputs required by the theory are taken from the Brownian dynamics simulations, there 
are no free parameters in this comparison. 

Section II presents the general theory for describing the nonequilibrium dynamics, while the BD simulation methods 
and united atom alkane potentials are described in the subsequent two sections. Theory and simulations are compared 
in Section V with no adjustable parameters. The final section discusses the potential application of the theory to 
more complex systems, such as proteins. 



II. NONEQUILIBRIUM AVERAGES FROM THE THEORY 

We consider a polymer consisting of N beads. Denote the 3N Cartesian coordinates of the polymer beads by 
r(i) = {ri(t),r 2 (t), . . . ,r 3N (t)} = {xi(t),yi(t), z\{t), x 2 (t), y 2 (t), z 2 (t), . . . ,xjd£Ly N (t), z N (t)}. The dynamics of this 
model polymer are assumed to be governed by the Smoluchowski equation,EJjt3 

aP(r ' <|r °' ^^(r)P(r,t|r o ,0), (1) 
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where P(r,t | i"o,0) denotes the probability that the coordinates of the N beads at time t are r, given that the 
coordinates of the N beads at time are ro, ks is Boltzmann's constant, T is the absolute temperature, U(r) is 
the potential energy, and £ is the bead friction coefficient. The only explicit degrees of freedom in these equations 
are the positions r of the polymer beads. The solvent affects the motion of the polymer through implicit frictional 
and stochastic forces. The frictional force is taken to be large enough that the motion lies is in the overdamped 
regime — that is, inertial effects may be neglected. We also neglect hydrodynamic interactions between different 
beads, but the general theory below may readily be applied including hydrodynamic interactions in equation (|^). 
The hydrodynamic interactions are omitted .here to reduce the computational labor but have been treated previously 
within the generalized Rouse-Zimm theory.li-3 The simplifying nature of these assumptions is useful for specifying 
a model for which it is possible to make rigorous, parameter-free comparisons between theory and "exact" results 
from computer simulations of equation (jl]). Moreover, this particular model contains many of the important physical 
features without being overly complicated. After understanding the application of the theory to this relatively simple 
system, the theory can be generalized to treat more realistic and complex systems. For example, the equilibrium 
theory has recentlv,been applied to describe the dynamics in alkane melts where the dynamics is governed by the full 
Liouville equations 

We are interested in describing the relaxation of a polymer from nonequilibrium starting conditions to equilibrium. 
In the example illustrated below, the nonequilibrium initial condition is the all-trans conformation of pentadecane, 
while in the final state, the pentadecane chain presents a large distribution of conformations. Denote by (f(t)) none q 
the expected value of some quantity f(t) for a system that starts from nonequilibrium starting conditions at time 
t = 0. The nonequilibrium average (/(t)) n oncq can be expressed in terms of probability distributions as 

(/(0>nonc q = J dr P(r , 0) J dr/(r)P(r, t | r , 0), (3) 

where P(ro, 0) denotes the probability that the coordinates of the N polymer beads at time are ro, P(r, t \ Tq, 0) is 
the formal solution to equation (Fy), 
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P(r, t | r , 0) = e tz? Wp(r, | r , 0) = e tv ^S(r - r ), (4) 
and S(r — ro) is the Dirac delta function. Thus, equation (|^) may be rewritten as 

(/W)no„c q = / dr o P(r o ,0) J drf(r)e tv ^6(r-r ) 

dr o P(r o ,0) J drS(v-v )e tc ^f(r) 

dr o P(ro,0)e t£ ( ro )/(ro), (5) 



where C is the adjoint of V, 
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The formal expression in equation (^) is rendered computationally tractable by introducing the eigenfunctions \& n 
and eigenvalues A„ of C, 

CV n = -\ n * n , 71=1,2,.... (7) 

Because e~ u ^ kBT ^ C is self-adjoint, the natural inner product is defined by 

/ \ _ J drexp[-U(r)/(k B T)]xiXj ^ 
Ui ' Xjj " JdreM-U(r)/(k B T)} W 

and coincides with the definition of the equilibrium average, 

, v v _ /<frexp[-t/(r)/(fc B r)] XtXj 

Jdrexp[-[/(r)/(fc B T)] ' W 

With this choice of inner product, a complete, orthonormal set {^n} of eigenfunctions can, in principle, be found. 
Thus, the /(ro) of equation (||) can be expanded in terms of the eigenfunctions {^n} as 

/(r ) - * n )*„(r ) = £</¥ B >* n (r ). (10) 

n n 

Substituting equation ( p"o| ) into equation (|5[) produces 

(/W>„o„ oq = / droP(ro,0)e t£ ( ro )^(/*„>* n (r o )-^e- tA "(/vI/ n ) / dr P(r , 0)* n (r o ) (11) 

n n 

Since the eigenvalue equation (|^) cannot be solved analytically for any realistic potential, the solutions are approx- 
imated using a basis set expansion, 

M 

*n«£C in &. (12) 

The coefficients Ci n are therefore determined by solving the matrix eigenvalue problem,0 

FC = SCA, (13) 

subject to the normalization constraint, 

C T SC = I, (14) 

where I is the unit matrix, 

Sij = (<t>i<j>j), (15) 
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and A is the diagonal matrix of (approximate) eig envalues X[, A' 2 , . . . , X' M . Substituting the Ci n and A^ obtained from 
the matrix eigenvalue problem into equation ( |ll[ ) produces the final equation for (f(t)) noncq as 

M / M \ / m r \ 

(/(*)>noneq ~ E (E ^«</&> J ( J2 C ^ J * P(r , 0)^ (r ) J e - A »*. (17) 

The nonequilibrium average in equation (|l7j ) is evaluated using as input only equilibrium information ((ftfri), F, and 
S) at the final temperature T and similar averages (J droP(ro, 0)</>j(ro)) for the initial condition. In our case, the 
initial condition is the all-trans conformation of pentadecane, and "averaging" over the initial condition simply involves 
evaluating cj>j for the all-trans conformation. However, equation (O) can be used for more general nonequilibrium 
starting conditions, such as an initial temperature Ti n i t (or pressurejdifferent from the final temperature T, in which 
case J droP(ro,0)(j>j(ro) is calculated by averaging cf>j over the Boltzmann distribution corresponding to Tj n i t . In 
fact, the initial condition could even correspond to a state (equilibrium or constrained) of the polymer in a different 
solvent. 

III. BROWNIAN DYNAMICS SIMULATIONS 

The BD computer simulations are more conveniently pursued by rewriting the Smoluchowski equation (Q) as a set 
of equivalent Langevin equations, 

where X*(t) is a Gaussian random variable with the properties, 

(X*(t)) = 0, 

(X?(t)r i (t / ))=0, fori>i', (19) 

(X*{t)X*{t'))=2^5 ij 6(t-t'), 

where Sij is the Kronecker delta function and 8{t — t') is the Dirac del ta. . fu nction. The Brownian dynamics simulations 
are performed using the algorithm of van Gunsteren and Berendsen,r°H 21 l 

n (n+l)=rdn) + ^At+±^(At) 2 +Kn, i = l,...,3N, (20) 

where At is the time step, Fi(n) = —dU(r) /dri is the force acting on the ith coordinate, F^(n) = [Fi(n) — Fi(n—1)]/ At 
is the finite difference approximation to the time derivative of the force, and \ in is a random positional displacement 
taken from a Gaussian distribution with zero mean and variance 2(fcsT/C)At. 

The calculation of the nonequilibrium averages (f(t)) none q from the BD simulations uses a scries of trajectories 
run at a temperature of 300 K starting with pentadecane in the all-trans conformation where the molecule lies in 
the x — y plane with the long axis parallel to the x axis and with the first three beads having the coordinates 
(0.00A,0.87A,0.00A), (1.26A,0.00A,0.00A), and (2.52A, 0.87A, 0.00A). Because all of the trajectories differ greatly 
from one another, obtaining good averages requires many trajectories, and the simulations use 9000 trajectories, each 
of 800 ps duration. 

The theory of section III calculates the equilibrium averages (• • •) from a single 400 ns trajectory run at a temperature 
of 300 K. The system is allowed to equilibrate for 20 ns, after which data are collected for 380 ns. 
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IV. POLYMER MODEL 



Pcntadecane is modeled as a string of 15i-GH2 or CH3 united atom groups. The potential U is, except for one 
modification, that of the GROMOS package£3 This potential has the general form 

U = [/bond + Wangle + ^dihedral + [/nonbond- (21) 

The bond length potential is 

14 

2 



14 1 

U hond = J2Mh-lof, (22) 



i=l 



where k is the length of bond i, lo — 1-53 A is the GROMOS equilibrium bond length for alkanes, and Kb is the bond 

length force constant. We use the force constant Kb = 160 kcal mol - A , which is smaller than Kb in GROMOS 
by a factor of 5. This change does not significantly affect the long-time dynamics but allows using a larger Brownian 
dynamics time step (At = 5 fs). The bond angle potential is 
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Y^-Ke^-do) 2 , (23) 



Wangle 2 
i=l 

where 0i is the ith bond angle, 9q — 111° is the GROMOS equilibrium bond angle for alkanes, and Kg = 
110 kcal moP 1 rad -2 is the bond angle force constant. The dihedral angle potential is 

12 

[^dihedral = ^ K$(\ + COS 30;), (24) 
1=1 

where fa is the zth dihedral angle, and K$ — 1.4 kcal/mol. The nonbonded potential [/ n0 nbond uses a 6-12 Lennard- 
Jones potential to describe the interaction between two nonbonded united atoms. More details are given in Ref. [2^. 

The solvent is modeled as a structureless continuum that exerts a viscous damping force on the pentadecane chain. 
The friction coefficient is calculated using Stokes' law £ = 67777 R, where the viscosity 77 is taken to be 1 cp and the 
hydrodynamic radius of a CH 2 or CH 3 united atom group is taken to be 1.5 A. 



V. RESULTS 

The only approximation in the theory is associated with the use of a finite basis set in equation (II4). Thus, we begin 
by studying the effectiveness of various basis sets for estimating the time evolution of the mean-square end-to-end 
distance R^ nd = ( x i5~ x i) 2 + (yi5~ yi) 2 + ( z i5~ z i) 2 - In contrast to the previously studied equilibrium time correlation 
functions of dynamical quantities that are odd in the position coordinates,EriJ the mean-square end-to-end distance is 
an even function of the position coordinates. Therefore, symmetry considerations dictate the use of different (even- 
power, as opposed to odd-power) basis functions in the mode-coupling theory calculation of (R 2 nd (t)) noneq . iTjhe scalar 
nature of R^ nd (t) implies that a general mode coupling basis may be constructed from the basis functions^ 

1, 



li • lj [V (distinct) pairs {i,j},i,j = 1, . . . ,14], 
(li • lj)(lk • lm) [V (distinct) pairs of (distinct) pairs {{i,j}, {k, m}}, i,j,k,m= 1, . . . , 14], 
(li • lj)(lk • lm)(l n • 1 P ) [V (distinct) triples of (distinct) pairs, 
{{h j}, {k, m}, {n,p}}, ij, k, m, n,p = 1, . . . , 14], 
etc., 
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where li denotes the ith pentadecane bond vector with components U x = JCj+i —Xi, h y = Ui+i — yi, and U z = Zi+i — Zi, 
for i = 1, . . . , 14. The ordering of pairs (and triples, etc.) is unimportant; for example, the pair {1, 2} is equivalent 
to {2, 1} since li • la = la • li, and hence only one of them is included in the basis set. The constant basis function 1, 



which produces a vanishing eigenvalue in equation (13), accounts for the non-zero value of (i?e nd (i))noncq at infinite 
times. 

The mode coupling theory is first applied using the relatively simple basis set, 

{1, li • lj [V (distinct) pairs {i,j},i, j = 1, • • • , 14]}, (Basis set I). 

Fig. 1 demonstrates that basis set I provides a reasonable estimate for the time evolution of (R 2 nd (t)) ncm cq whose 
numerically "exact" dynamics is represented by the BD simulation. A more accurate approximation of (-Rend(*))noneq 
requires expanding basis set I to include more basis functions of the appropriate symmetry, such as the tetralinear 
products (li • lj)(lk • lm)- However, the total number of tetralinear products scales as N 4 , where N is the number of 
monomers. Thus, inclusion of all tetralinear functions makes the basis set unmanageably large even for the relatively 
small "polymer" pentadecane. Clearly, all the tetralinear products cannot generally be retained in the basis set, 
nor are all these basis functions relevant to the long time dynamics. Those tetralinear products that contribute the 
most to the long time dynamics are determined using a selection technique similar to that introduced previously for 
computing equilibrium time correlation functions .EjTj 

Based on prior treatments of equilibrium time correlation functions, it proves convenient to work with basis functions 
constructed as linear combinations of individual bond vectors that provide a first approximation to the collective chain 
motions. An approximate set of eigcnfunctions, denoted as the first order GR (generalized Rouse) modes, of C may 



be obtained from equations (12) - (|lq ) with a basis set consisting of all linear functions in the bond vectors, 

{l lx [i = l,...,U],ky[i = l,...,U],kz[i = 1.....14]}. 

The first order GR modes are thus linear combinations of the bond vector components. Note that the linear basis 
functions and the corresponding first order GR modes do not contribute to the expression in equation (|T^) for 
(^cnd(O)noncq due to symmetry considerations. Products of first order GR modes, however, may be used as providing 
a first approximation to (-RcndW) noncq- Due to the invariance of C under interchange of the Cartesian directions 
x, y, and z, the first order GR modes need only be computed for one direction, say x, and the GR modes in the 
other two directions follow trivially by symmetry. Let ^iJ, ^2x, ■ ■ ■ , ^i2c be the first order GR modes obtained from 
solving equations (|l^) - (|l6|), and let their corresponding (approximate) eigenvalues , A 2 , . . . , A^ be ordered in 
increasing magnitude. Thus, the first eigenfunction provides the first order estimate to the GR eigenfunction 
with the slowest decay, called the slowest mode, provides the first order estimate to the GR eigenfunction of the 
next slowest decaying mode, and so on. Let ^%y, ■ > ■ • • > *i2, be the GR modes obtained by replacing all the x's 
in \f4* , i ■ ■ ■ j ^(n-i) x wrtn ^' s ' ana - *i» >*2* ' ■ ■ ■ j ^i4z be the corresponding GR modes in z. It is convenient 
to introduce the notation *>"^[i = 1, 2, . . . , 14] to represent the vector (^r ;', , )• We expect that the slowest 
decaying first order eigenfunctions contain the most information about the long-time dynamics. Thus, rather than 
including all possible tetralinear products in the basis set, we only include tetralinear .products constructed from the 
slowest decaying GR modes. This produces the basis set which contains the functionst.il 



{1, • [V (distinct) pairs {i,j},i,j = 1, . . . , 14] 



• ^^X*^ ■ *W) [V (distinct) pairs of (distinct) pairs {{i, j}, {k, m}}, i, j = 1, . . . ,Q}}, 



(Basis set Hq), 

where < Q < 14. The Q = case reduces basis set IIq = o to basis set I, while using Q = 14 is equivalent to including 
all possible tetralinear products in the basis set. Note that the size of basis set IIq scales roughly as N 2 + Q 4 . Hence, 
basis set IIq is significantly smaller than the full basis set containing all possible tetralinear products if Q <C N. 
Fig. 1 demonstrates that basis set IIq with Q = 4 produces an improvement in the long time dynamics predicted by 
the theory, though for short times below 100 ps, basis set I is in better agreement with the simulations. Increasing 
Q beyond 4 does not produce any significant improvement. Greater accuracy of the theory therefore requires the 
inclusion of hexalinear and higher order products in the basis set. However, the size of the basis set grows extremely 
rapidly as higher order functions are included, even when a procedure similar to that in basis set IIq is used. In order 
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to select those higher order that contribute the most to the long time dynamics, we now test an alternative sorting 
procedure first described in Ref. |J for treating the long time dynamics of small flexible peptides. 

The sorting procedure orders basis functions according to the sum of the (approximate) eigenvalues of the GR modes 

in the tetralinear and higher order (even) products of GR modes. For example, A,- 1 "* + Xj + A^ + A™ represents the 
first order relaxation rate for the tetralinear product (\&- ■ ■ V&m*)- Since we expect that the functions 

most relevant to the long time dynamics are composed of the slowest-decaying GR modes (smallest approximate 
eigenvalues), the basis set retains the products with the smallest eigenvalue sums. While first order eigenvalues may 
not provide an optimal estimate for the relative importance of different products of GR modes, ptejaous experience 
suggests that this method is sufficient for obtaining excellent equilibrium time correlation functionsJj'u The first order 
eigenvalue sorting procedure produces the basis set which is represented as 

{1, [V (distinct) pairs {i,j},i,j = l,...,U}, 

higher order GR mode products (up to 14 products) with the R smallest eigenvalue sums}, 

(Basis set III fi ). 

Fig. 1 indicates that basis set III^ with R = 500 produces further improvement in the theory. Increasing R beyond 500 
results in numerical instabilities in .the calculation, though previous experience indicates that using longer equilibrium 
simulations rectifies the problem.™ 

The discrepancy in the "infinite" time R^. nd in Fig. 1 is due to the difference between the "infinite" time (i.e., 
equilibrium) (R^. nd (t)} noncq from the 9000 nonequilibrium BD simulations and the (i? 2 nd ) from the 400 ns equilibrium 
BD simulation. This difference arises from statistical fluctuations in calculating averages, and the magnitude of these 
fluctuations can be estimated by a simple calculation. The variation in i? 2 nd is comparable to the size of i? 2 nd itself 
- about 2 ran 2 . In a BD simulation, the "memory" of previous values of i? 2 nd is about 100 ps (see Fig. 1), so a 
400 ns BD simulation yields about 4000 independent samples of i? 2 nd . Because i? 2 nd is not normally distributed, the 
central limit theorem is not strictly applicable.^ Nevertheless, the theorem can provide an estimate for the order of 
magnitude of the fluctuations in calculating the average (i? 2 nd ) as 2 nm 2 /v4000 ~ 0.03 nm 2 . This value is consistent 
with the discrepancy in Fig. 1 and is also consistent with the variation in the equilibrium average (-R 2 nd ) between two 
different 400 ns BD simulations (data not shown). 

Fig. 2 demonstrates that the relaxation of -R 2 nd to equilibrium is highly nonexponential at long times. This is in 
contrast to equilibrium dipole time correlation functions, where the exponential decay at long times is predicted very 
accurately by the mode coupling theory.Q The nonexponential behavior of (i? 2 nd (i))noneq at long times may account for 
the difficulty in obtaining highly accurate calculations from the theory. Nevertheless, it is encouraging that the mode 
coupling theory gives reasonable estimates for (R^. nd (t)} noncq , and, furthermore, correctly predicts that the relaxation 
at long times is nonexponential. 

VI. DISCUSSION 

The computational effort required to calculate (i? 2 nd (t)) n oneq directly from Brownian dynamics simulations is enor- 
mous and requires orders of magnitude more computer time than is necessary for simulating equilibrium time cor- 
relation functions. The absence of time translational symmetry during a nonequilibrium relaxation precludes the 
computation of this process from a single Brownian dynamics trajectory. Instead, an ensemble average must be per- 
formed over a huge number of nonequilibrium trajectories. For example, 9000 trajectories are necessary to obtain with 
sufficiently small noise the relaxation to equilibrium of (-R 2 nd (i)) n oneq for pentadecane from the all-trans conformation. 
The accumulation of these nonequilibrium trajectories requires more than one month on an SGI Power Challenge. 
Moreover, this simulation is far simpler than that needed for most nonequilibrium processes. The direct simulation 
of a T-jump process would necessitate averaging over the distribution of both initial and final conditions and would 
require at least an order of magnitude more than the 9000 trajectories used for the simulation with constrained initial 
conditions. Clearly, any realistic study of the unfolding of a protein by computer simulations is currently unfeasible 
and will remain so even in the more distant future. In contrast, the theory requires as inputs only equilibrium aver- 
ages which may be obtained readily from a few Brownian dynamics trajectories. Furthermore, the theory can easily 
handle many different types of nonequilibrium processes. Treating the relaxation of pentadecane from the all-trans 
conformation requires only equilibrium averages for the final state. To describe a T-jump process using the theory, 
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equilibrium averages must be computed for both the initial state and the final state — merely a doubling of the 
computational effort. 

The successful description by the theory of the "unfolding" of pentadecane indicates that the theory proposed 
here may provide a significant contribution to the study of protein folding and unfolding. Understanding how a 
globular protein folds requires detailed characterization of partially organized intermediates formed during the folding 
process. Information about such intermediates is becoming available through many complementary experimental 
studies involving stable partially folded states, protection from hydrogen exchange, disulfideJjormation in oxidative 
refolding, effects of amino acid substitutions, and peptide analogues of folding intermediates.E£l Nonetheless, no clear 
unifying view of the nature of folding mechanisms has yet emerged, and many uncertainties over fundamental issues 
remain. 

One reason that protein folding remains an unsolved problem is that the experimental techniques used to study 
this process do not provide structural information with atomic detail about the multitude of folding pathways and 
intermediates. Instead the folding process is generally characterized experimentally by more global properties, such as 
the degree of compactness of the protein, the percentage of secondary structure compared to the native state, different 
hydrogen exchange rates of various residues, etc. The real time observation of protein folding with, for example, NMR, 
circular dichroism, or fluorescence techniques, also provides information about coarse-grained features of the folding 
intermediates and pathways with limited spatial resolution. Consequently, the folding process is frequently analyzed 
in terms of the free energy as a function of a heuristic one-dimensional reaction coordinate which is not directly 
derived from microscopic motions and thus only reflects global qualitative features of the folding process. Computer 
simulations provide limited information as well because only short time scales are accessible to realistic simulations. 
Furthermore, the large amount of computation required by such simulations allows the study of only one (or at best 
a few) trajectories from which it is difficult to extract general characteristics of the folding process. 

The theory presented here potentially lacks many of the limitations of experiments and computer simulations. 
In principle, the theory can predict the time evolution of any interatomic distance in a polypeptide and can thus 
specify the average structure of the polypeptide at any time during the folding or unfolding process. However, the 
practical application of the theory to protein folding may emerge only after further testing and development. For 
example, presently it is unclear how the accuracy of the theoretical predictions scales with system size. Therefore, it 
is desirable to test the theory for systems larger than pentadecane. A natural choice would be the helix-coil transition 
in a simple polyalanine helix since the computational intensity of nonequilibrium simulations for a small helix is not 
significantly greater than that for pentadecane. Thus, nonequilibrium BD simulations are feasible for the unfolding 
of the small helical peptide, thereby again enabling no-parameter tests of the theory before applying the theory to 
even larger systems for which the nonequilibrium simulations are prohibitive. The recent treatment for the more 
complicated equilibrium dynamics of small peptides with the mode coupling theory lends support to the belief that 
the nonequilibrium generalization presented here is adequate for the task of describing the unfolding o£-a small protein. 

A possible direct experimental test of the theory is provided by the work of Perkins and coworkers.E3 They measure 
the relaxation of a single DNA molecule from a highly stretched nonequilibrium state to equilibrium using optical 
tweezers to hold the DNA molecule in place and fluorescence labeling to follow the motion of the DNA molecule. The 
relaxation of the end-to-end vector is highly non-exponential, with a fast initial decay followed by a slower decay. The 
simple Rouse and Rouse-Zimm models do not fully explain the data; this is not surprising since modeling DNA as a 
string of beads connected by springs is clearly an oversimplification. The mode coupling theory provides a framework 
for treating DNA more realistically. Much work remains to be done, however, before a direct comparison can be 
made between the experimental data of Perkins and coworkers and the predictions of the mode coupling theory. The 
extension of the theory to include the solvent would be an important first step, and the treatment of hydrodynamic 
interactions is possible following Ref. |f4] Also, the applicability of the theory to highly charged polymers, such as 
DNA, must be tested. Nevertheless, the predictive power of the mode coupling theory for describing the dynamics of 
a variety of polymers suggests that such generalizations would be fruitful. 
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FIGURE CAPTIONS 



Fig. 1. The relaxation of the mean-square end-to-end distance (R^. nd (t)) nonC q of pentadecane to equilibrium from the 
all-trans conformation. The solid line represents the "exact" result calculated from the average of 9000 nonequilibrium 
BD simulations, while the various dashed lines represent the results of mode coupling theory calculations using different 
basis sets, as described in the text. 

Fig. 2. Same as Fig. 1, except that the (Rg nd (t)) noneq axis has a logarithmic scale. 



9 




100 200 300 400 500 600 700 800 

Time (ps) 



Figure 1 Tang; J. Chem. Phys. 
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Figure 2 Tang; J. Chem. Phys. 
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